Nonlinear propagation of light in structured media: 
generalized unidirectional pulse propagation equations 
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Unidirectional pulse propagation equations [UPPE, Phys. Rev. E 70, 036604 (2004)] have pro- 
vided a theoretical underpinning for computer-aided investigations into dynamics of high-power 
ultrashort laser pulses and have been successfully utilized for almost a decade. Unfortunately, they 
are restricted to applications in bulk media or, with additional approximations, to simple waveg- 
uide geometries in which only a few guided modes can approximate the propagating waveform. 
The purpose of this work is to generalize the directional pulse propagation equations to structures 
characterized by strong refractive index differences and material interfaces. We also outline a nu- 
merical solution framework that draws on the combination of the bulk-media UPPE method with 
single- frequency beam-propagation techniques. 



I. INTRODUCTION 

Computer simulations in the field of nonlinear op- 
tics have been playing an important role in understand- 
ing ever more extreme regimes in light-matter interac- 
tions. Dynamics of ultrashort, high-power laser pulses 
is one particular field, which motivated significant ef- 
fort and concomitant progress in numerical methods de- 
signed for optics at femtosecond time scales. One can 
say that optical filamentation played the role of a cata- 
lyst for the development of a number of pulse propaga- 
tion models, which made detailed studies of extremely 
nonlinear regimes a possibility. However, the most accu- 
rate pulse propagation models remain restricted to bulk 
media, both gaseous and condensed, while waveguiding 
structures have to be treated with more significant ap- 
proximations. 

The purpose of this article is to put forward a theoret- 
ical framework, which will allow the implementation of 
simulators capable of handling pulse propagation regimes 
characterized by the following four attributes: 

A) Structures with strong refractive index contrasts. 

B) Directional long-distance wave propagation. 

C) Strong waveform reshaping, both in time and space. 

D) Extreme spectral dynamics, with resulting spectra 
often encompassing more than an octave in fre- 
quency. 

This combination is rather difficult to handle numerically. 
For example, there exists a wealth of work (e.g., 
utilizing the beam propagation method (BPM), which 
are designed for regimes A and B, and can incorporate 
certain weak nonlinearities @, but are restricted to nar- 
row spectral regimes. Time-domain beam propagation 
methods have been developed (e.g. 0,11]), though they 
concentrate mainly on linear regimes. Direct Maxwell's 
equations solvers [Il-Qjj are well-suited to regimes A, C, 
and D, but are prohibitively expensive for simulating 



long-distance pulse propagation. Simulators based on 
the Unidirectional Pulse Propagation Equation (UPPE) 
fu| and other types of one-way propagation equa- 
tions can cope with attributes B, C, and D, but 
require additional approximations to simulate waveguid- 
ing structures, such as hollow-slab leaky waveguides 24 \ 
or nonlinear nano- waveguides (25l - [27| . In other words, 
methods suitable for the combination A+B+C+D have 
yet to be developed. 

In this paper, we present a step in this direction and 
describe a generalization of the UPPE, which can be ap- 
plied to nonlinear structured media with strong differ- 
ences between refractive indices of the constituent mate- 
rials. Departing from the wave equations, we derive an 
auxiliary evolution system. This is used to find projec- 
tion operators that extract forward and backward prop- 
agating components of the field from an arbitrary optical 
field waveform. These operators transform the auxiliary 
system into a coupled forward-backward pulse evolution 
system that is exact and accounts for structured media. 

Key to the approach we present is that the pulse evolu- 
tion equations are cast in a form which makes it possible 
to combine proven numerical methods. More specifically, 
nonlinear interactions can be treated by ordinary differ- 
ential equation (ODE) libraries, the same way it has been 
done with UPPE-based simulators [28[ ; the linear propa- 
gator can be treated by tapping the rich knowledge base 
of BPM, and in particular the techniques developed for 
wide-angle BPM (WA-BPM) (see Refs. [H, for early 
formulation, and Refs. [3ll . H2| for examples of various 
Pade approximated propagators). 

The remainder of this paper is organized as follows. 
First, in Sec. HH we give a brief summary of the UPPE 
model. In Sec. IIII1 this model is generalized and a cou- 
pled forward-backward pulse evolution system derived. A 
unidirectional propagation approximation is applied and 
the resulting equation transformed into a form analogous 
to a bulk-medium UPPE. Considering a homogeneous 
medium finds the generalized equation reduces back to 
UPPE. Section IIVI outlines a strategy for a numerical 
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solution of the generalized system. We summarize and 
discuss future directions in Sec. [V] 



II. UNIDIRECTIONAL PULSE PROPAGATION 
EQUATIONS: SUMMARY 

The main purpose of this paper is to generalize the 
UPPE framework. For the reader's convenience and 
ease of reference, we briefly summarize the corresponding 
equations. 

In a homogeneous medium characterized by a dielectric 
permittivity e(w), a pair of coupled UPPEs are exact 
[T3I [H[ and can be written in the form 



d z E + (k±,0J,z) — +ik z E + (k±,uj, z) + 
,2 



(1) 



E 

8=1,2 



IbJ 



2enc 2 fc 



■P(kj_, U), Z) 



2e c 2 h 



-J(k±,oj,z) 



d z E^(k±,uj, z) = —ik z E^(k±,u>, z) — 



(2) 



-P(k±,u>,z) - 



2eoc 2 fc ; 



■ J(fc_L, OJ, z) 



These equations describe the evolution of E±(k±, uj, z), 
which are the spectral (Fourier) representation of the 
electric field, k^ = {k x ,k y ,0} are the transverse wave 
numbers and k z is the z— component of the wave- vector 



k = |fc x , k y , k z = ^Ju 2 e(uj)/c 2 



1,2 



K V I 



(3) 



which satisfies the dispersion relation k 2 = u> 2 n 2 (ui)/c 2 . 
The two polarization vectors e s {k±,uj) are orthogonal 
to k and to each other, but otherwise can be chosen 
freely. The superscript _L denotes the transverse part 
(i.e., x,y) of the corresponding vector. Equations (jTJ) and 
([2]) are mutually coupled through the nonlinear medium 
polarization P(k±,uj, z) and current density J(fcj_, w, z). 
These responses are functionals of the electric field. They 
are normally specified in the real space and time repre- 
sentation: 

P(r ± ,t) = P({E(x,y,t)}) , J(r x ,t) = J({E(x,y,t)}) . 

It has to be emphasized that the system of Eqs. (JXJ) and 
@ is exact and together with the V • D equation (which 
can be used to obtain the z component of the field if 
needed), is equivalent to Maxwell's equations. However, 
as with direct Maxwell's equations solvers, it would be 
difficult to solve in its entirety, i.e., including forward and 
backward propagating waves. In practice, the unidirec- 
tional propagation approximation is assumed, and the 
medium response is calculated solely from the forward 
propagating waveform: 



Under this approximation the system reduces to a sin- 
gle UPPE, Eq. (TT]). For details of the numerical so- 
lution, the reader is referred to Ref. |28j]. Here we 
only point out that the native representation suitable for 
numerical implementation relies on spectral amplitudes 
j4 S) _|_(fcj_, uj, z), which only change with z due to nonlinear 
interactions with the medium. They are related to the 
electric field through the linear propagator e lk '( k ^< k v u ') z : 



E^(k x ,u,z)= J2 et^,+ (fc±,^K Mfc ^ )z . (5) 

8=1,2 



The corresponding UPPE equation, 

i,ip~ lk ^ z 

d z A s , + (k^,u},z) = 



2e c 2 k z 



(6) 



e s ■ [iuP(k±,u,E + {z)) - J(k±,u,E+(z))]. 



P{E) 1 J{E)^ P(E + ),J(E^ 



(4) 



constitutes a large system of ODEs. This is the repre- 
sentation in which it is solved numerically. Because the 
medium response is calculated in the real time represen- 
tation at each spatial point, spectral transforms in both 
directions have to be invoked multiple times when the 
right hand side of the ODE system is evaluated. 

The main limitation of the UPPE approach is that 
it is restricted to homogeneous media. Weakly guiding 
structures can be included as part of the polarization 
response, but geometries with strong material contrasts 
and interfaces cannot be efficiently simulated. Reference 
[l3j shows derivation of the UPPE system for waveguid- 
ing structures, but its implementation requires knowl- 
edge of the full system of electromagnetic modal fields, 
which is impractical to obtain even for geometries that 
admit exact solutions. Therefore, waveguiding scenar- 
ios can only be simulated under additional assumptions, 
which require that the fields can be described as a super- 
position of a few guided or leaky modes of the structure, 
whatever field configurations evolve. 

The practical limitation of the UPPE, and in fact of 
all other uni-directional pulse propagation methods orig- 
inates in the identification of the forward and backward 
propagating waves resting on the usage of a reference ho- 
mogeneous medium. To elucidate this, let us consider the 
following example. Let e(w) represent a chosen homoge- 
neous background, and let x(ri> w ) be such that e+x( r ±) 
gives the actual permittivity of the structure. The vari- 
able part x( r J-) can be treated within the polarization 
term of Eq. ([6]), P(r) — eQx{ r ±-)E+(r±_). Now note 
that even if x is constant throughout space, the propaga- 
tion constant of a plane wave predicted by Eq. i.e., 
k z (uj,k±) + x^ 2 /(2c 2 k z (uj, k±)), is only a second order 
Taylor approximation to the exact plane-wave propaga- 
tion constant ^ui 2 [e(uj) + x(w)]/c 2 — kj_. On the other 
hand, it is straightforward to show that if we retained 
both coupled UPPEs ([]]) and ©, the resulting propaga- 
tion constant would be exact independently of the choice 
of the reference e(u>). We thus see that if x( r J_) varies 
in space, no matter how we select the reference medium 
e(w), even the linear problem is not solved exactly by a 
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single UPPE. This means that if we wish to use the uni- 
directional approximation for nonlinear interactions, we 
must find a way to marry it with an exact propagation 
description in the linear limit. This in turn implies that 
no part of the refractive index variation in space should 
be included in the polarization term of the propagation 
equation. 



III. DIRECTIONAL PULSE PROPAGATION 
EQUATIONS 

In this section, we generalize unidirectional pulse prop- 
agation equations to situations with material interfaces 
parallel to the propagation direction z and with strong 
refractive index differences between materials that com- 
prise the structure. Central to this task will be the 
ability to extract the true forward and backward prop- 
agating components of the total electromagnetic field. 
The main deviation from the method described above 
is that the pulse propagator native representation will be 
mixed; we will retain the spectral representation of the 
frequency /time dimension, but will use the real space 
representation for the transverse spatial dimensions x, y. 



A. Model of a nonlinear, structured medium 

Consider a non-magnetic, isotropic dispersive medium, 
with the dielectric permittivity e(x,y,uj), which only de- 
pends on coordinates x and y and angular frequency uj. 
We assume there are no free charges or currents. The 
constitutive relation for all media will be written in a 
form using polarization to account for all properties ex- 
cept the linear e(w) 

P = -p(x,y,{E(x,y,t)}) . (7) 

We assume that an algorithm is given that computes po- 
larization from a given history of the electric field vector 
E(x, y,t) at a specified point [x,y]. The first two argu- 
ments of P are meant to indicate that this algorithm can 
depend, through the medium properties, on the trans- 
verse location [x, y] but not on the longitudinal coordi- 
nate z. The concrete functional form of P is unimpor- 
tant for the present purposes, but for a specific example, 
the reader can think of the instantaneous optical Kerr 
effect in which the local index of refraction changes pro- 
portionally to the square of the electric field vector. As 
the medium is isotropic, the polarization direction follows 
that of E: 

PKerr(x,y, {E(x,y,t)}) 

= 2e n 2 (x,y)(E 2 x + E 2 + E 2 z )E. (8) 

Here, nzix, y) stands for the nonlinear index, which as it 
indicates, may depend on location. Other models of light- 
matter interactions that have been used in simulations 



are described in articles on filamentation [3J,l35| and Ref. 
|28j shows methods for their numerical implementation. 

To keep notation simple, we will not use current den- 
sity explicitly. In general, nonlinear interactions with the 
medium can be equivalently formulated cither in the po- 
larization or current density language, so this means no 
loss of generality. In numerical simulations, using both 
current and polarization may actually be convenient, and 
it only requires trivial extension of our results. 

B. Fields in terms of analytic signals 

In numerical simulations, it is often easier to work with 
analytic signals of the electric field. Here we use analytic 
signals to represent all real quantities (E(t),P(t)). For 
example, the electric field is obtained as a real part of its 
analytic signal: 

E = Re{E(x,y,z,t)}, (9) 
which has its spectrum restricted to positive frequencies: 

E(x,y,z,t)= / dujE(x,y,z,oj)e- iut . (10) 
Jo 

Here, and in what follows, we will distinguish be- 
tween temporal and spectral representations of functions 
through their respective arguments t and uj. Because 
the only time we need the representation of the electric 
field in the time domain is when we compute the non- 
linear medium response (i.e., polarization), we will work 
mostly in the spectral representation. 

C. Derivation of Directional Pulse Propagation 
Equations 

Our departure point is the wave equation for the elec- 
tric field, accompanied by a constraint in the form of the 
divergence equation. 

VV-J5- V 2 E= (eE+^-Pj , V-D = (11) 

The divergence equation deserves a note. While we have 
assumed no free charges and currents, high intensities 
can lead to medium ionization and subsequently to elec- 
trons drifting away from their parent ions. However, our 
treatment aims to describe femtosecond pulses. On such 
a short time scale we can safely assume that even when 
ionization occurs, the positive and negative charges do 
not have enough time to separate, and the average lo- 
cal charge remains zero. Therefore, using the divergence 
equation, 

V-D = e eV-E + e E- Ve + V • P = 0, (12) 

V • E can be expressed in terms of the nonlinear polariza- 
tion divergence and transverse electric field components 
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as follows 



-V-E = -S±- V±e+— V-P (13) 
e e e 

The transverse (x, y) part of the wave equation is thus 
rewritten to separate the linear and nonlinear terms 



d zz E ± =LE ± +N ± [E}. 



(14) 



Here, the linear operator L is related to the corresponding 
Helmholtz equation (for a fixed angular frequency). It 
acts only on the transverse electric field vector 



ur 1 
LE ± = —e(r ± ,w)E ± +A ± E ± +W-E ± .W ± e. (15) 
c z e 



The nonlinear operator N acts on E±, but in general, 
also depends on the E z component 

uj 2 1 
N[E} = — =P(E)+V — V-P(E). (16) 
e c z e e 

We will address how to obtain E z later. For now, let us 
assume that it can be calculated once E x<y are known. To 
obtain propagation equations for E±(z,x,y,u), we first 
introduce auxiliary field amplitudes, effectively doubling 
the number of variables used to describe the electric field: 

Ei(z,x,y,u) = E+(z,x,y,uj) + E^(z,x,y,w), (17) 

with 



E+ = A+(z,x,y,Lu)e + ^ 
E~ = 4-(z,a;,y,w)e- lCz , 



(18) 



where i = x,y and £ stands for a parameter to be chosen 
freely. Clearly, since no £ appears in Maxwell's nor the 
wave equations, physical observables must not depend on 
the concrete choice of £, and this "gauge invariance" will 
become manifest when we arrive at our final result. Wc 
will term £ a reference wavenumber to emphasize the fact 
that it has no physical meaning by itself. 

It is also important to keep in mind that the positive 
and negative wavenumber parts E i of the field are, in 
general, not the forward and backward propagating por- 
tions of the total waveform. So far we have not restricted 
how fast A i can change with z. In principle, they could 
evolve so fast that their variation would completely over- 
ride the exponential factors e ±II ' z accompanying them. 
That is why both E i can contribute to waves propagat- 
ing in the positive and negative z direction (see Ref. [36[ 
for how this occurs). 

By representing a single function Ei{z,x,y,ui) as a 
combination of two functions Ef (z, x,y,u), we have 
added artificial degrees of freedom. These will be taken 
back by requiring that Ef satisfy a relation of our choice. 
Concretely, we impose an additional constraint in the 
form 

e+^ z d z A+{z, x, y, u) + e~ iiz d z Ar ( z , x, y,u) = (19) 



The rationale behind this constraint is exactly the same 
as in the variation of constants method for differential 
equations. Namely, this representation eliminates the 
second derivatives when one evaluates d zz E. Because 
of the constraint, the first derivative simplifies to 

d z E, = t((E+ - E^) (20) 

and the second derivative to 

d zz E t = -C 2 E t + i(e +l < z d z A+ - Ke~< z d z A- (21) 

Using this in the wave equation (fTH) together with the 
constraint of Eq. (TT9l , we obtain the evolution equations 
for the auxiliary amplitudes A^ 



d z Ai 



±i 
2C' 



(L - ( 2 )E ± + N[E] 



(22) 



To evaluate the right hand side of this system, E± is ex- 
pressed in terms of Af, and E z is subsequently obtained 
from the z-component of the wave equation. Using Eq. 
(|2U)l . the latter can be written as follows 



Kd x (E+ -E-) + tCd y (E+ - E-) - ^P Z {E) 



d X xE z + d yy E z H —E z 



(23) 



If not for nonlincarity, this is an inhomogeneous 
Helmholtz equation that determines E z in terms of the 
transverse field components. The part of the polariza- 
tion component P z which is nonlinear in E z , is usually 
very small and therefore normally neglected. Should one 
not be satisfied with such an approximation, the above 
equation can be solved by iteration. For example, in Kerr 
media for intensities typical of femtosecond filaments, a 
single iteration already gives an accurate result. 

The next step consists in identifying the parts of the 
electric field waveform which propagate in the positive 
and negative directions along the z axis. The resulting 
equations become more intuitive when expressed in terms 
of auxiliary iJ-fields: 

d z E+ = +i(E+ + JL [(i _ f)E ± + N[E]\ (24) 



d z E~ 



2C 



(L-C 2 )E ± +N[E] (25) 



This system, completed by Eq. (f2"3"]l. is equivalent to the 
wave equation with the divergence constraint, and can be 
solved in principle. However, in this form, it poses two 
problems. First, in general, it would require very short 
propagation steps in order to resolve both the forward 
and backward propagating waves. Second, the physical 
input conditions for simulations are normally given such 
that the problem to solve is a boundary value problem 
rather than an initial value problem. The latter point be- 
comes evident when we realize that it is only the forward- 
propagating field component that is specified at z = 
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(e.g., at the laser output). The second condition is that 
the backward propagating field is zero at z — > oo (i.e., 
at the far end of a laboratory). Such a boundary value 
problem would be rather difficult to solve. Fortunately, 
in many cases the backward propagating wave can be ne- 
glected. We shall therefore derive the beam propagation 
equations that account for such a situation, but in the 
process shall identify the forward and backward propa- 
gating fields (and see that they are, in general, different 
from i? ± ). 

In matrix notation, the propagation equations read 



El 




E- 



Having separated the linear and nonlinear part of the 
evolution operator, we are in the position to determine 
the forward and backward propagating parts of the field. 
This division will be defined with respect to the linear 
system. In the spirit similar to Ref. 37], two projector 
operators can be constructed from the Helmholtz opera- 
tor L and its square root L^: 



V b 



4C 
4C 



+ (C + ^) 2 +(L-C) 

+(i-C 2 ) +(C + ^) 2 



(27) 
(28) 



It is straightforward to show that these operators have 
the expected properties of projectors, in particular they 
are idempotent 

V% =V F V%= V b (29) 

and they constitute a unity decomposition 

V F + V B = 1 , V f Vb = V f Vb = 0. (30) 

These projectors also commute with the linear evolu- 
tion operator in Eq. (|2^|) . and direct calculation shows 
that their eigenvectors have propagation constants cor- 
responding to forward and backward modes propagating 
in the linear system. Thus, we can use these projectors 
to obtain the true forward and backward propagating 
field components. If the total field is given in terms of 
the auxiliary amplitudes E^, then the forward portion 
of the wave is obtained as 



E F = (1 l)V F 



El 
E7 



(E] 



■e: 



- L-k(El - E-)\ . (31) 

This expression contains the reference wavenumber, 
which might suggest that it depends on our artificial split 



of the field in Eq. However, because of Eq. (j20l) . 

the above expression reduces to 



Ep 



1 



:E> 



-L-*d z E ± 



(32) 



2 2 

and the backward propagating component is obtained as 
1 



E r 



:E A 



\L-id z E A 



(33) 



As they must be, these forward/backward amplitudes are 
independent of the reference wavenumber. Our aim is to 
express the pulse evolution equations in term of these 
directional fields. Because the projector operators are z- 
independent, the simplest way is to apply them directly 
to the propagation equations. In other words, we need 
to compute 



d z E\ 



d z E\ 



= (1 1)V F 
= (1 l)V B 



d z E+ 
d z E- 

d z E+ 
d z E7 



(34) 
(35) 



Inserting the right hand side from Eq. (|26[) and using the 
projector properties of Eqs. (|2U)) and (13TH) . we obtain a 
pair of coupled equations for the forward and backward 
fields 



d z E F = +iVLE F 
d z E B = —i"sft,E B 



2v L 

i 



2VL 



N±[E F + E B ] 
Nj_[E F + E B ] (36) 



This is a generalization of the coupled pair of Unidirec- 
tional Pulse Propagation Equations ((T|) and ((2]). We 
show later that in a homogeneous medium, for which 
we have an explicit expression for the square root of the 
Helmholtz operator, these equations reduce to UPPEs as 
they should. 

Similar to the case for bulk media, Eq. ([3"6l is not 
the best for numerical implementation. We transform 
this system into a form analogous to that of bulk UPPE 
© such that we can adopt the same numerical solution 
strategy. Toward this purpose, we use amplitudes which 
will exhibit evolution only if some nonlinearity is present: 

El = e+^AKz) , Ef = e- l ^ z Af(z) (37) 
In this representation, Eq. (1361) reads 

Q zA F = J^ e -iVt*ft ± [ e +iVL* A F + e -V£* A Bl 



2v L 



d z Af = 



+iV Lz jq _ L [ e +»V Lz j^F _|_ -iy Lz a B 



2v L 



A 1 



(38) 



This shows explicitly that the forward and backward 
propagating waves are mutually coupled in the nonlin- 
ear terms. It is obvious that for strong nonlinearity, our 
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forward-backward projection loses its intended meaning, 
because it can renormalize and couple waves propagating 
in both main directions [38]. Thus, we arrive at a point 
where we must adopt an approximation which will allow 
us to reduce the full system to a single unidirectional 
equation. 



D. Unidirectional propagation approximation 

Our final step is to adopt the unidirectional approxi- 
mation, where we assume that the nature and strength of 
nonlinearity is such that only negligible backward prop- 
agating fields are generated. Then, the nonlinear term 
can be approximated as 



N ± \e +iVLz A F + e -^ Lz A B ^ 



N±[e 



(39) 



and the system can be restricted to only the forward- 
propagating field: 



d z A F (r±,u,z) = +- 



— iV L;. 



Z A F ] 



2v L 



(40) 



This is the sought-after generalization of the Unidirec- 
tional Pulse Propagation Equation. As expected, the 
structure of this system is completely analogous to the 
bulk UPPEs of Eq. ^ , with the exception that the linear 
propagator is formally expressed in terms of a Hclmholtz 
square root operator, instead of plane-wave expansion. 
The most pronounced difference is that Eq. (|40|) is na- 
tively represented in the mixed representation. It retains 
the spectral treatment of the time dimension and with 
that, it preserves the ability to treat chromatic and non- 
linear properties of the material exactly. On the other 
hand, the transverse dimensions are represented in real 
space, which is the natural choice for the implementation 
of the linear propagator in a structured medium with 
strong refractive index variations. 



then use Eq. (JST 



8 Z E F (k ± ,u,z) = +ik z E F (k ± ,Lo,z) + —N ± [E F ] , (42) 



and express N in terms of polarization: 

■ 2 1 



d z E l 



ik 7 E* 



— yP(E) kk ■ P(E) 



(43) 

Only the transverse components in this equation consti- 
tute the evolution system, but in this full-vector form, it 
is easy to see that the operator acting on the polariza- 
tion term produces the transverse part of the nonlinear 
response, namely 



uj z 1 

— ttP(E) kk ■ P(E) 

e c z e e 



kk- 



P(E). 



The projector operator in the square brackets can be re- 
placed by a sum over projectors on the polarization vec- 
tors e, 



1 - 



kk- 



E 



(45) 



Using this in Eq. (|44|) and inserting it into Eq. (j43|) . we 
obtain 



d z E F (z,uj, k±) = ik z E F - 



2e c 2 fc ; 



^eie s -P(E), (46) 



which is identical to the homogeneous medium UPPE of 
Eq. ([1]) (with the current density term omitted). Thus, 
as it must, the generalized pulse propagation equation 
([4"U| passes this sanity check and reduces to the UPPE if 
the medium is homogeneous. 



E. Special case: reduction to UPPE in a 
homogeneous medium 

Before going into how this pulse evolution equation can 
be solved numerically, let us illustrate how it reduces to 
the well-known bulk UPPE for a homogeneous medium. 
First, we recall that for a homogeneous medium, we know 
that plane waves are eigenfunctions of the Hclmholtz op- 
erator L and that in the plane-wave representation, the 
linear propagator reduces to multiplication by a phase 
factor given by the propagation constant k z (u, fcj_): 

^ — i^f~Lz ^ — ik z (uj,k±)z 

It is therefore sufficient to Fourier-transform Eq. (|40"]) 
from the (x, y) space to the transverse wavenumber space 
(k Xl k y ) to obtain 

d z Al{k u uj,z) = +^- e -^ z N ± [e + ^ z A F l (41) 



IV. NUMERICAL SOLUTION STRATEGY 

In this section we sketch, in broad strokes, an approach 
for the numerical solution. It builds on the ODE-based 
method for solving UPPE systems and combines it with 
a wide-angle beam-propagation solver used to evaluate 
the linear propagator exp(iL 1 / 2 z). 

The core of the pulse propagator of Eq. (|40l) is an ODE 
system, with z being the independent variable. The equa- 
tion is evaluated at every transverse spatial location r± 
and frequency w while being incremented along the prop- 
agation direction z. During a single ODE step, the right 
hand side of Eq. (|40p has to be evaluated multiple times 
at different values of z which are subject to the choice 
of the specific ODE algorithm. Because the integration 
is normally executed with an adaptive integration step, 
one cannot determine beforehand at what specific z lo- 
cations the term [exp(iL 1 / 2 z)A] needs to be computed - 
an algorithm is needed to evaluate the right hand side 
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for any small value of z. For the ODE solver, we use the 
open source Gnu Scientific Library (GSL), but any im- 
plementation with the following capabilities can be cho- 
sen. One necessary feature of a suitable ODE library is a 
driver for adaptive step control, with a robust algorithm 
monitoring the accuracy of the numerical solution. An- 
other necessity is that the library contains methods which 
do not require Jacobian evaluation, because such meth- 
ods are not suitable for UPPE-like ODE systems [28| . 
We typically employ the Runge-Kutta-Fehlberg method, 
however, another useful ODE library feature is the capa- 
bility to switch easily between different integration meth- 
ods. Regarding its structure and method of solution as an 
ODE system, the generalized propagation equation does 
not differ from an ordinary UPPE. What is different is 
the implementation of the linear propagator. 

Because the linear propagator is diagonal in angular 
frequency, this task is equivalent to a set of uncoupled 
beam-propagation problems. In other words, the action 
of exp(zi 1 / 2 z)'!/; only requires one independent BPM-likc 
update for each uj resolved in the simulation. This por- 
tion of the algorithm is therefore "embarrassingly paral- 
lel," with perfect balance and no inter-dependencies be- 
tween calculations performed for different angular fre- 
quencies. There are many wide-angle BPM methods 
available, and any of them can be utilized, in principle. 
For instance, one can evaluate the linear propagator by 
a Pade approximant. Defining (3 2 (lu) = ui 2 e(ui)/c 2 , the 
dominant part of the Hclmholtz operator, one writes 



LAz 



n 



X + a k 
X + b k 



(47) 



The coefficients a k , b k depend on Az and are chosen as 
to reproduce the Taylor expansion of the left hand side. 
For example, 



Ai + (i - (3Az)X 
Ai + (i + (3Az)X 



ipAz 



(48) 



is second-order accurate in X with the error scaling as 
~ /3AzX 3 . Various higher order approximations can be 
constructed in the same spirit. If the operator X acts in 
a non-trivial way along both spatial dimensions x, y, it 
is often further split into "one-dimensional" components 
so that the resulting matrices are band-diagonal. 

Similar techniques can be used to compute the inverse 
square root of L that acts on the nonlinear response term 
in Eq. (|40[) . However, this operator can be approximated 
by L 1 / 2 w u>n(ui)/c as is usually done in filamentation 
simulations [281 ]. This is sufficient unless the spatial pro- 
file of the nonlinear polarization becomes "focused" to 
wavelength scale. 

With the linear propagator implemented as a "BPM- 
based plug-in," the solution proceeds in steps with two 
stages: 

1. Call the ODE solver. One integration step is exe- 
cuted that updates the current A(r±,u>, z) into the 



new A(r±,uj,z + Az). The ODE solver algorithm 
invokes computation of the right hand side of Eq. 



'LSz 



2/3 



LSz 



which contains two applications of the linear prop- 
agator [e.g., Eq. ([i5])] for a sub-step ±5z. Be- 
hind the scenes, the solver determines the maxi- 
mum step Az possible on a global scale, since some 
parts of the grid may contain finer features, and re- 
quire shorter integration steps than others. Unlike 
the fully spectral UPPE, the length of the integra- 
tion step the ODE solver is permitted to take is 
bounded from above by the maximum step allowed 
by the BPM method used for the linear propagator. 

2. Re-align the spectral amplitudes. The point along 
z at which A and E amplitudes coincide can, of 
course, be chosen arbitrarily. It is advantageous to 
renew this synchronization point after each ODE 
step such that A and E coincide at the beginning 
of the ODE step. This is achieved by 



A, <- 



'LAz A F 



(49) 



which amounts to yet another application of the 
linear-problem propagator to the current solution. 
Naturally, Az must be obtained from the ODE 
solver as the actual length of the last adaptive in- 
tegration step. This repeated re-alignment step is 
normally implemented in bulk-media UPPE solvers 
as well, but there implementations without it are 
possible, in principle. Here, it is crucial that the 
step length in the linear propagator is kept small, 
and application of Eq. (|49| ensures that 5z is al- 
ways smaller than the maximal step allowed in the 
ODE solver. 

In a nutshell, the above procedure describes the stan- 
dard UPPE solution method modified in two ways: First, 
a BPM-based propagator is utilized for the (short-step) 
linear advancement of the optical field, and second, real- 
space representation of transverse dimensions is retained 
at all times. 



V. SUMMARY 

We have presented a generalization of the Unidirec- 
tional Pulse Propagation Equation suitable for structures 
characterized by material interfaces parallel to the pulse 
propagation direction and by strong differences between 
the properties of the constituent materials. While the 
main result of Eq. (|40l) is somewhat intuitive, we show a 
rigorous derivation based on identification of the forward 
and backward propagating wave components. These are 
expressed in terms of projection operators [Eq. (I28[) ] akin 
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to those we have previously used in bulk media 12] . They 
allow expression of the generalized UPPE in terms of the 
linear propagator, and they "isolate" the nonlinear inter- 
actions with the medium, such that the evolution is de- 
scribed in terms of spectral amplitudes which only evolve 
due to non-zero nonlinearity. 

The generalized UPPE uses a mixed representation: 
spectral for the time/frequency dimension and real-space 
for the transverse (to the direction of propagation) di- 
mensions. The linear propagator can be based on one of 
the many available beam-propagation methods. The con- 
crete choice of the method will depend on the given ge- 
ometry. For example, the important case of nearly plane- 
parallel wave guide structures [H, H3] can be treated by 
a one-dimensional WA-BPM combined with the plane- 
wave expansion in the free-propagating direction. In- 
dependently of the chosen BPM approach, the numeri- 



cal solution strategy developed previously for bulk-media 
UPPEs can be used with relatively minor modifications. 

There is an increasing interest in extreme nonlinear op- 
tics confined to wave-guiding structures of different kinds. 
It is therefore expected that our results will find appli- 
cation in various implementations of efficient pulse prop- 
agation solvers, especially situations in which both the 
geometry of the structure and waveform reshaping due 
to nonlinear interactions play important roles. 
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